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Abstract 

Mean field models (MFMs) of cortical tissue incorporate salient, average fea- 
tures of neural masses in order to model activity at the population level, 
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thereby linking microscopic physiology to macroscopic observations, e.g., 
with the electroencephalogram (EEG). One of the common aspects of MFM 
descriptions is the presence of a high dimensional parameter space capturing 
neurobiological attributes deemed relevant to the brain dynamics of interest. 

We study the physiological parameter space of a MFM of electrocortical 
activity and discover robust correlations between physiological attributes of 
the model cortex and its dynamical features. These correlations are revealed 
by the study of bifurcation plots, which show that the model responses to 
changes in inhibition belong to two archetypal categories or "families" . After 
investigating and characterizing them in depth, we discuss their essential 
differences in terms of four important aspects: power responses with respect 
to the modeled action of anesthetics, reaction to exogenous stimuli such as 
thalamic input, distribution of model parameters and oscillatory repertoires 
when inhibition is enhanced. Furthermore, while the complexity of sustained 
periodic orbits differs significantly between families, we are able to show how 
metamorphoses between the families can be brought about by exogenous 
stimuli. 

We here unveil links between measurable physiological attributes of the 
brain and dynamical patterns that are not accessible by linear methods. 
They instead emerge when the nonlinear structure of parameter space is 
partitioned according to bifurcation responses. This partitioning cannot be 
achieved by the investigation of only a small number of parameter sets and 
is instead the result of an automated bifurcation analysis of a representative 
sample of 73,454 physiologically admissible parameter sets. Our approach 
generalizes straightforwardly and is well suited to probing the dynamics of 
other models with large and complex parameter spaces. 

Keywords: bifurcation analysis, brain dynamics, mean field model, 
electroencephalogram (EEG), anesthesia, thalamus 



1. Introduction 

Understanding brain function is one of the longstanding and unresolved 
quests of science. Despite the efforts and contributions of a wide range of 
scientists from different disciplines, we are still very far from a comprehen- 
sive account of what the brain does, and how and why it does it. The 
formidable hardness of this problem is embodied in the complexity of the 
structure and organization of this organ, which supports interactions over a 
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wide range of temporal and spatial scales. One of the scientific approaches 
that aims to shed some light on the complexity of brain activity is the use 
of realistic mathematical models and computer simulations: constrained by 
known physiology and anatomy, these theories try to recreate and understand 
the patterns emerging in the electrical activity of the cortex. In particular, 
so-called mean field models (MFMs) [l7|, 33, 22] try to capture the results 



of neuronal interactions through spatial and/or temporal averaging, which 
then describe the population activity of regions of cortex in a parsimonious 
manner. Accompanied by a sufficiently flexible parametrization, MFMs take 
into account existing uncertainties in cortical anatomy and physiology in an 
effective way. 

Despite the apparent biological simplicity of MFMs, such theories produce 
physiologically and behaviorally plausible dynamics. Since their birth forty 
years ago, they have been employed successfully to investigate a large number 
of cognitively and clinically relevant phenomena, such as the acti on p harma- 



cological agents have on the activity of cortical tissue |42|, |43|, |25|, [14J, |26|, |45 
the origin of evoked potentials 
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rhythm |37|, |29 



ship to cognitio n |44j. 131 



seizures 



41, 32 



the sleep-wake cycle and circadian 
the emergence of gamma band oscillations and its relation- 
and the onset and characteristics of epileptic 
c.f. Ref. lal for a recent review. We use here 
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Liley's MFM [22| , which is capable of reproducing the main spectral features 
of spontaneous (i.e., not stimulus-locked) EEG, in particular the ubiquitous 



alpha rhythm [20j. Firstly, this model is biologically constituted: all state 
variables and parameters can be constrained on the basis of existing anatom- 
ical and physiological measurements in the literature. Secondly, it supports a 
rich repertoire of behaviors both physiologically relevant and dynamically in- 
teresting. For example, parametrically widespread, robust chaotic activity of 
various origins has been found 0,139, 13], and multistability, i.e., the presence 
of various coexisting dynamical regimes, has been demonstrated. Multista- 
bility has been speculated to correspond neurobiologically to the formation 
of memories 0] . 

in this paper we link nonlin- 



In the spirit of the dynamical approach 15 



ear electrical activity and neurobiologically significant attributes of cortex. 
To this end, we consider a representative sample of parameter sets that have 
previously been found to generate physiologically plausible behavior Q . This 
sample contains 73,454 sets and for each of these we computed the bifurcation 
plots when varying two parameters related to inhibition, i.e., the qualitative 
changes of recurrent activity patterns, when neuronal inhibition is altered. 
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It turns out that we can sort the sets into two distinct "families" of dynam- 
ical behavior. These families are found to correlate with EEG signal power 
and responses to anesthetics, whereas family membership is determined by 
specific neurobiological parameters. 

The paper is divided into three main sections. First, we briefly intro- 
duce Liley's ordinary differential Equations (ODEs). We next discuss the 
theoretical and numerical tools employed in some detail, in particular how 
one arrives at a systematic bifurcation analysis procedure to show the qual- 
itative changes of the model solutions when inhibition is varied. Then, we 
present the main distinctive features of the families, such as their responses 
to the simulated induction of anesthetics and their correlations with model 
parameters of interest. We are also able to show how exogenous agents, 
i.e., input from the thalamus to the cortex, can induce dramatic changes 
in those patterns and stimulate transitions from one type of family to the 
other. That, in particular, provides a compel ling example for the modulation 



thalamus is thought to exert on the cortex |35j. All these relations cannot 
be discovered with standard linear or nonlinear analyses of the physiological 
parameter space, and represent the main result of this work. A discussion of 
open problems concludes the paper. 



2. Neuronal mean field equations 

Liley's MFM aims to provide a mathematically and physiologically parsi- 
monious description of average neuronal activity in a human cortex, with spa- 
tially coarse-grained but temporally precise dynamics. One excitatory and 
one inhibitory neuronal population, respectively, is considered per macrocol- 
umn, which is a barrel-shaped region of approximately 0.5 — 3 mm diameter 
comprising the whole thickness of cortex (thus ~ 3 — 4 mm deep). Corti- 
cal activity is locally described by the mean soma membrane potentials of 
the excitatory (h e ) and inhibitory (hi) neuronal populations, along with four 
mean synaptic inputs J ee , Jj e , I e i, and la. These inputs convey the reciprocal 
interaction between neuronal populations, where double subscripts indicate 
first source then target (each either excitatory e or inhibitory i). The con- 
nection with measurements is through h e , which is linearly related to the 
EEG signal j28[. Lumped neuron populations are modeled as passive RC 
compartments, into which all synaptically induced ionic currents terminate. 
According to population types (j, k) = e, i, synaptic activity drives the mean 
soma membrane potentials from their resting values. The equations for h e 
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and hi are given by 



dh- h eq -h(t) h eq -h(t) 

where h r e and hi are mean resting potentials, and r e and r, are the mem- 
brane time constants of the respective neuronal populations. The reversal 
potentials of the transmembrane ionic fluxes mediating excitation and inhi- 
bition are given by h e e g k and h^, respectively Note that the synaptic inputs 
are weighted with +1 (excitatory l e ^) and —1 (inhibitory 1^) at the resting 
potential of the respective excitatory or inhibitory neuronal population, and 
that these weights then vary linearly with voltage. 

The mean synaptic inputs describe the postsynaptic activation of ionotropic 
neurotransmitter receptors by presynaptic action potentials, arising from the 
collective activity of neurons both nearby and distant. The time course of 
such activity, based on well-established experimental data [38], is modeled by 
a critically damped oscillator driven by the mean rate of incoming excitatory 
or inhibitory axonal pulses. We thus have, for k — e,i: 



\ 2 

+ lek) hk{t) = F ek J ek e jiV^Se [h e (t)) +p ek (t) + e jfc(t)j , (3) 

+ Hk\ hk{t) = r ife7ifc e|Arf fc »S i [/i i (t)]+p ifc (t)} , (4) 



d 
~dt 
d_ 
~dt 

where the terms in curly brackets correspond to sources of the axonal pulses 
from three origins: local, i.e., in the same macrocolumn of the cortex N^Si, 
arriving through long-range, excitatory cortico-cortical connections from other 
macrocolumns e fc, and extracortical, i.e., primarily of thalamic origin pi k . 
For subsequent simplicity we assume the absence of any extracortical in- 
hibitory input, i.e. pn~ = 0. N^ k quantifies the strength of anatomical pop- 
ulation connectivity. The maximal postsynaptic potential (PSP) amplitude 
To, occurs in the target population k = e, % at time after the arrival of 
the presynaptic spike from the source population I = e,i. A schematic illus- 
tration of the architecture of interactions in the Liley model can be found in 

Fig. m 
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Figure 1: Architecture of Liley's mean field model. Two separate model macro- 
columns are shown, each containing one excitatory and one inhibitory neuronal popula- 
tion. Note that long-range connections are exclusively excitatory and that self-couplings 
correspond to connections of neurons of the same type within the local populations. 

Local mean soma potentials are nonlinearly transformed to mean neu- 
ronal population firing rates with a sigmoidal function 

S k [h k (t)] =S r r x {l + exp 

where \i k and <j k indicate the firing thresholds and their associated standard 
deviations for the respective neural population. The chosen form is a com- 
putationally convenient approximation to the error function, which results 
when one assumes a Gaussian distribution of firing thresholds. 

The final part of the model concerns the propagation of excitatory axonal 
activity in long-range fibers. In the full formulation, the following inhomo- 
geneous, two-dimensional damped wave equation is used: 

<P ek (x,t)=N: k v 2 A 2 S e [h e (x,t)] . (6) 

Two core assumptions form the basis of this equation: first, action potentials 
traverse much larger distances in the cortico-cortical fibres than the local 



-V2 



hk{t) - fik 



<7k 



(5) 
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(intra-columnar) axonal fibres, and thus their conduction delays can no longer 
be considered negligible. Second, cortical areas further away share fewer 
long range connections. Given the average number of long range excitatory 
connections onto a local neuron N™ k , experiments suggest that the fraction 
thereof that stems from distant neurons decays approximately exponentially 
with distance, with a scale constant A. Activity is then modeled as spreading 
omnidirectionally with constant conduction speed v. We have also made 
the additional assumption that long-distance propagation dynamics is not 
differentiated according to excitatory and inhibitory targets, i.e., A ee = A e j = 
A and v ee = v e i = v. 

In the following, we force the term V 2 e fc — > and thus consider only 
bulk dynamics for our model cortex. Furthermore, we neglect the spatially 
inhomogeneous effects of long-range fibre systems. Thus, Eq. ([6]) now be- 
comes 

7 \ 2 

( +oo) <f )ek (t) = N: k u J 2 S e [h e (t)} , (7) 



dt 

with bj = vA. We see that Eq. ([7]) has turned into an inhomogeneous critically 
damped oscillator. This is in contrast to simulating an unconnected local 
population only, which would involve dropping Eq. fl6]) altogether. However, 
it is clear that the steady state solution to Eq. ([7]), (p e k(t) = N" k S e [h e {t)\, 
simply changes Nj? k — > N" k + Nj? k in Eq. Thus, the steady state of bulk 
oscillations corresponds to local activity with increased connectivity. 

The fourteen ODEs of Eqs. (j2])- (SJ) and (JTj) are studied in the present pa- 
per. We note for the interested reader that a thorough discussion of realistic 
axonal velocity distributions in MFMs has been recently published by Bojak 
and Liley j^] and that Liley et al. 22] contains more in-depth discussions 



of the physiological and mathematical arguments that lead to Eqs. ([2])- ([6]). 
Also, parameters used in the equations and their physiological ranges are 
shown in Tab. [TJ 



3. Methods 

Parameter sets and bifurcation parameters 

A large number of parameter sets is considered in this study, represent- 
ing a significant sample of the physiologically admissible parameter space of 
Liley's MFM. To generate such realistic parameter sets, the properties of the 
model have been studied around fixed points of linearized version of Eqs. ([2])- 
dHD- Biologically plausible selection criteria were then invoked concerning the 
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properties of the simulated EEG and the firing rates of their associated neu- 
ronal populations. Limits on the admissible ranges of parameter values were 
also applied, so that they are consistent with well-established neurophys- 
iological features of mammalian cortex. Details on the methods and the 
constraints applied in this search can be found in Bojak and Liley j^j. As 
a result, a collection of 73,454 physiologically meaningful parameter sets is 
available and considered in this study as representative of physiologically re- 
alistic parametric instantiations of model cortical activity in Liley's MFM. 
Bifurcation theory is an essential tool for discussing the qualitative changes 
that solutions of dynamical systems undergo under variations of their param- 
eters. Loosely speaking, this type of analysis provides a picture similar to a 
phase diagram for physical systems, so that changes in the properties of the 
solutions (i.e., bifurcations) may be considered as a generalization of phase 
transitions [36]. This allows us to explore parametric boundaries between 
qualitatively different activities. 

One of the a priori strongest motivations for the use of this approach is 
represented by the failure of standard methods of statistical linear analy- 
sis, which are unable to supply relevant information about the model's 29- 
dimensional parameter space. An example of this is illustrated in Fig. [2j As 
a result of a Principal Component Analysis (PCA) over the parameters of 
the 73,454 sets, the percentage of total variance explained by the first ten 
PCs is plotted, individually and cumulatively. Essentially, if the parameter 
space presented strongly linear correlations among the variables, it would 
be possible to identify a limited number of components capturing the vast 
majority of the variance, indicating that the significant degrees of freedom 
in the space were in fact much less than its dimension. This is clearly not 
the case, since the first ten components cumulatively account for only about 
42% of the total variance, with low and very similar singular contributions. 
It is not possible to find a linear combination of parameters that expresses 
the "bulk" properties of the parameter space. If we want to extract valu- 
able statistical information, a different, non standard method of analysis is 
needed. On this basis, the aim of this paper is to highlight correlates between 
dynamical "snapshots" of the model, i.e. two parameters bifurcation plots, 
and spectral and parametric properties of Liley's MFM. 

In choosing the parameters to vary for our bifurcation analysis, we are 
guided by the existence of literature pointing to inhibition being a sensitive 
locus of control of brain dynamical responses. For instance, it is well-known 
that pharmacological agents alter cortical activity primarily through act- 
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ing on inhibition |34|. Factors that modify inhibitory PSP amplitudes and 
anatomical coupling strengths of the inhibitory neural population are thus 
selected as continuation parameters. This is in line with previous work j3] 
and implies changes in Equation H] for the inhibitory PSP amplitudes 

Tie — > R^ie ) Tjj — > RTu , (8) 

and local inhibitory-inhibitory connectivity 

N&-+kNi, (9) 

where R and k are chosen as our bifurcation parameters. 

This choice is not exclusive and other parameters or factors could have 
reasonably been varied, possibly highlighting other poignant features that 
have not emerged through continuations in R and k. Nonetheless, these two 
parameters allow biologically plausible control: R changes the amplitude of 
the inhibitory PSPs for both the excitatory and inhibitory populations, which 
for example general anesthetics and sedative agents can affect to a large de- 
gree. At the same time inhibitory-inhibitory connection strengths are altered 
via k, which changes the impact on the inhibitory feedback loop largely re- 
sponsible for generating specific oscillations j22|. Note that since pik = in 
Eq. m the effect of inhibition on the excitatory population scales precisely 
with R, but on the inhibitory population with R ■ k. Thus, for example, 
setting k = 1/R would simulate an exclusive impact on the excitatory popu- 
lation. Essentially then, we can control overall inhibitory "potency" with R 
and "specificity" with k. 

One and two parameters bifurcation analysis procedure 

There are two types of bifurcation plots that will be discussed: one pa- 
rameter (lpar) plots in R and two parameters (2 par) graphs in R and k. Let 
us first explain how the diagrams are obtained: we integrate Eqs. and 
(J7j) for a chosen parameter set up to an equilibrium solution for R = 1 and 
k — 1, i.e., in the absence of changes to inhibition. Initial conditions are 
given at t = by h e (0) = /^(O) = — 70 mV, while all the other state vari- 
ables are set to zero. Once an equilibrium solution is found, it is continued 
in R and a lpar plot is obtained. In all the i?-plots there are always two 
types of codimension one (codl) bifurcations present: saddle-node and Hopf 
points, respectively. The presence of Hopf bifurcations gives us some infor- 
mation about the character of sustained oscillations produced by the action 
of inhibition, which we will discuss below. 
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The saddle-node and Hopf bifurcations are subsequently continued in 
parameters R and k. In our diagrams, lines describing the loci of saddle- 
nodes are labeled sn, and similarly for the Hopf points hb. 

Physiological limits for minimal and maximal parameters intervals are 
set at R m in = k min = 0.75, R max = 2.00 and k max = 1.25. Higher or 
lower values of R, k are considered to alter inhibitory PSP amplitudes and 
inhibitory-inhibitory connectivity strengths either to unphysiological values 
or to unusually large or small figures, which, for example, do not capture the 
commonly observed effects of pharmacological agents on the cortex. However, 
as discussed below, continuations are performed till particular features of 
the bifurcation topology become apparent. Typically, this occurs far outside 
of the physiological range. Bifurcation analysis is performed here with the 



widely used continuation software packages AUTO-07P and MATCONT 10, 

8- 



4. Results 

Two parameters bifurcation patterns 

To establish our findings, we proceeded with a two-step analysis. First, 
2par bifurcation plots for a reduced sample of sets were obtained and ana- 
lyzed, visually and numerically. Analysis of plots in R, k for 405 randomly 
selected sets out of the 73,454 in the generated sample revealed the existence 
of two unique, and clearly distinct, categories of patterns which recurrently 
appear. Our primary aim in this subsection is to describe these archetypal 
patterns and discuss their shared and distinct features. In the next subsec- 
tion, a third parameter (jp e i) is varied and this shows that the two families are 
"transformable": it is possible to metamorphose one family into the other, 
by varying exogenous stimuli. Once the two families are characterized, gen- 
eralization through the whole batch is straightforward and global features of 
the physiological parameter space naturally emerge. 

The general character and topology of the 2par plots for Liley's MFM 
is reproduced in Fig. El These represent two archetypal _R/c-plots for so- 
labelled Family 1 (Fl) and Family 2 (F2) types of sets. The major difference 
among the two families is represented by how the sn lines are organized and 
divide the parameter space in areas with different dynamical properties. In 
general, because of the form of the MFM's equations, a region bound by sn 
lines has always three equilibria, whereas the region outside has only one. 
For Fl, see Fig. sn lines travel almost parallel for the whole range of 
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the space, and do not form cusp points (c), dividing the space essentially in 
three major regions: one inside sn lines and two outside. In contrast, F2, see 
Fig. [3)3, is characterized by the presence of two cusp points Ci and C2, so that 
the region containing three equilibria is the union of two separated "wedge- 
shaped" areas, with cusps as their vertices. The way hb and sn lines interact 
is, for the majority of cases, by so-called fold-Hopf (fh) points, which appear 
as tangencies between such lines. Other cases we saw but are not present in 
Fig. |3]happen when hb lines connect two Bogdanov-Takens (bt) points. These 
instead correspond in the diagram to points where a hb line terminates on a 
sn line. In general, these bt points occur on the left, lower wedge of the sn 
lines, close to fhi for Fl and F2. We will give examples of these bifurcations 
in the next subsection. 

For both families, there are also other cases that slightly differ from Figure 
|3l although the structure of the sn lines never changes. Variations are minimal 
and local, i.e. one of the fh points can be missing or a small difference in the 
shape of the hb lines can be present, but the essential structure of the plot is 
not altered. What instead shows high variation, independently of the family, 
is the extension of the plots, i.e., the values of R and k at which bifurcation 
points appear. It is important to realize that the physiological region is 
typically a tiny area of the plots, as indicated by the rectangles in Fig. |3j 
There are plots, as those shown in Fig. [3j where the whole pattern is contained 
within 3 — 4 orders of magnitude of the bifurcation parameters, but there 
are also a large number of patterns for which hb and sn lines extend up to 
immense values, like R ~ 10 9 . This is an interesting, important point: scales 
vary widely, but the nature and structure of the patterns are invariant, and 
this is appreciated only if plots are considered in their complete extension in 
Rk-sp&ce. Moreover, it appears that similar bifurcations maintain the same 
type of normal forms among equivalent patterns, no matter what scale is 
involved. This implies that the dynamics in the proximity of bifurcations 
is topologically equivalent, i.e. qualitatively very similar, among plots in 
the same family. Thus, parameter sets in the same family share not only 
the topology of the saddle-node bifurcation lines, but also yield qualitatively 
similar dynamics close to the codimension two points. 

Thalamic input and pattern family metamorphosis 

The magnitude of extracortical (thalamic) input to inhibitory neurons 
p e i is chosen as a third continuation parameter for investigating the transi- 
tion between families. Our choice of p e i has an intuitive meaning: unlike T ie 
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and Ta, whose variation represent endogenously driven effects on the cortex, 
changes in the thalamic input capture the action of exogenous causes. Fur- 
thermore, p e i excites the inhibitory population and hence can be expected to 
affect significantly the dynamical repertoire of families based on variations 
of inhibitory PSPs. 

Let us discuss a prototypical transition from Fl to F2. To follow our 
description the reader will need to refer to Fig. H] . We start at a very large 
and unphysiological value of p e i = 26,000/s, indicated in Fig. HA, where a 
representative of Fl is present. Compared to the (standard) diagram depicted 
in Fig. [3j two fh points on the left wedge of line sni are not yet formed. 
When the input is decreased to p ei = 14, 000/s in Fig. HJ3, thus approaching 
the physiologically plausible range of 0/s < p ei < 10, 000/s, the hb branch 
extends towards lower values of R, creating the fhi point, as found in the 
standard Fl type in Fig. [3j Analogously, the fh 3 point is created at larger 
values of R and k. At the same time, on the right wedge sn 2 , two cusps Ci 
and c 2 appear in Fig. Hf3. This happens in a so-called swallow tail bifurcation 
(see, e.g., Ref. [lj]). All equilibria belong to the same branch, meaning that 
they are connected at the saddle-node bifurcations, at which the branch is 
folded. At this swallow tail bifurcation, a pair of saddle-node bifurcations 
appears on the unique branch of equilibria. Notice in the inset that the fh 2 
point is now on the right wedge of the little triangle formed by the two cusps 
Ci and c 2 . Moreover, the single branch hb has lost its previous elliptical 
shape, and has produced intersections. However, these are not bifurcation 
points since only one unique branch of hb lines is present, and the diagram 
still represents a type of plot belonging to Fl. The emergence of the two 
cusps is not considered a major change, since the topological characteristic 
that we choose in discriminating Fl from F2 is the global shape of the sn 
lines. When separate sn x and sn 2 do not intersect but keep almost parallel, 
as in Fig. HA and B, a diagram of Fl type is still present. 

As the extracortical input is further decreased to p e i = 7, 000/s in Fig. HP, 
the triangle of line sri2 with the two cusps Ci and c 2 expands and moves 
to the left of the diagram, intersecting the opposite sr\ 1 branch (see also 
the inset in Fig. HP). Now c 2 is overlapping the line of saddle-nodes sn!, 
and the cusp Ci has crossed sni and is above it. Notice, in particular, the 
vicinity of the blue and green lines of sn in the inset: this represents a stage 
very close to the transition of Fl into F2. In fact, at smaller p e i values, 
the cusp c 2 further lowers and separation between the two branches of sn 
further diminishes, so that increasing portions of sn! and sn 2 become nearly 
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tangent. This condition persists up to a value of p e i where sn lines finally 
exchange, so that at p ei = 6, 800/s the two cusps no longer belong to the same 
branch, as Fig. HP and its inset show. The cusp Ci is now part of sni and C2 
belongs to S112: a second swallow tail bifurcation has taken place, signaling 
the appearance of an F2 type of diagram. In other words, two regions of 
the i?fc-space having triangular shapes have formed and overlapped, and C2 
now tends to move towards the right as p e i is decreased, separating further 
from Ci. A further decrease of the thalamic input will eventually result in 
the same plot as in Fig. [3)3, i.e., the prototypical F2 with two disjoint V and 
A-shaped curves of saddle-nodes. 

A physiologically relevant observation can be made at this point. One 
can investigate the effect of varying the excitatory thalamic input to either 
inhibitory {p e i) or excitatory (p ee ) cortical neurons. It turns out that there 
is a relatively robust family-specific effect of thalamic excitation, which is 
opposed for inhibitory and excitatory cortical targets: Fl (F2) is the pre- 
ferred family at high (low) values of p e i, whereas Fl (F2) is the preferred 
family at low (high) values of p ee . The incidence of family metamorphosis for 
variations within the physiological interval of thalamic input 0/s — 10, 000/s 
is quite strong: about 56% of the 73,454 sets undergo a transition when p ei 
is varied and approximately 42% when p ee is. Once again, inhibition seems 
to be somewhat more effective in controlling the global dynamical properties 
of the cortex, in this case as driven by excitatory exogenous inputs. 

Changes within families due to thalamic input 

To complete our analysis of 2par plots, let us make some remarks on the 
other instantiations of Fl and F2 that have not been illustrated yet, but 
do occur in the 405 analyzed sets. These diagrams are only slight, local 
variations of the ones shown so far, and still exhibit the global common 
characteristics typical of the respective families. As an example, patterns 
belonging to Fl very similar to Fig. [4A are present, but they miss the fh.3 
point. This type of diagram is still considered as part of Fl, because the 
structure of its sn lines does not change. If we start looking at the transition 
between Fl and F2 from this particular plot, rather than from Fig. |4A, the 
fate of the pattern as p e i is decreased is the same, i.e., it metamorphoses into 
an F2 diagram. 

There is one final example of possible variations, regarding the appear- 
ance of Bogdanov-Takens (bt) points, which we anticipated in the previous 
subsection. Remarkably, these changes can be accessed as well by varying 
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the same parameter p ei : exogenous agents do not only produce transitions 
between families, but can also trigger modifications within the same family. 
We study here the annihilation of a couple of bt points, by increasing p e i for 
the same parameter set employed in Fig. HI Once again the reader is asked 
to look closely at Fig. [5] and its insets to follow our description. If we start at 
a low, physiological p ei = 3, 720/s, two Bogdanov-Takens points, bti and bt 2 , 
are present just below fhi and we want to demonstrate how they disappear 
for larger p e j. These bts occur on the left wedge of the sn line culminating 
in the cusp point C\. Notice also that the branch hbi is responsible for the 
point fh 2 close to Ci, in analogy with the F2 in Fig. |3j The similarity between 
the two plots is evident. As the thalamic input is increased to p ei = 3, 724/s 
in Fig. 03, a couple of generalized Hopf (gh) points appear on hbi. These 
points signal a change in the hb line from subcritical to supercritical, i.e. the 
Hop f bifurcation associated to the branch switches from a hard to a soft type 
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li Pei is made even larger, the two hb lines come closer and closer, up to a 
value of the input where they glue together and separate, as shown in Fig. \&P 
for p ei = 3, 729/s. After this, hbi acquires a large portion of what previously 
was hb2, with the latter now extending over a much smaller distance. The 
split occurs in between gh 1 and gh 2 , each of which now belongs to a separate 
Hopf branch: the tiny hb 2 reaches from bt! to bt 2 and the long hbi bends 
in proximity of gh 2 and continues to larger R and k values. Notice in the 
inset of Fig. that also hb 2 (in green) bends close to gh x and heads back to 
bti, where it stops. Then, when p e i is further incremented, see Fig. [5p with 
p e i = 3, 735/s, gh 2 annihilates and the smaller branch hb 2 further shrinks, 
with the bt points getting closer and closer as p e i grows. In the end, these two 
bt points coalesce, and we are left with only a single branch, as in Fig. [3B. 
We note that this line of hb delimited by bti and bt 2 only exists at positive 
values of R and k. For the interested reader, we point out that the plot 
in Fig. [5A is responsible for chaos in an unphysiological range of R and k, 
and its route through a so-called homoclinic doubling cascade is analyzed in 
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Let us conclude with some observations. First, the tiny hb 2 branch exists 
for a very limited interval of p e % and this is reflected in its very low occurrence 
in the 405 sets: it was found only twice. Second, an identical mechanism 
for Fl is also present, where the creation or annihilation of bt points takes 
place on the left srii lines in Fig. [3A. Finally, the diagrams that have been 
presented for the annihilation of bt points in F2, the similar plots for Fl (not 
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shown), and the plots in Figs. E] and S] with their small, local variations, are 
the only ones that have emerged in our analysis of the reduced batch of 405 
sets. Within the previously discussed assumptions, we conclude that these 
diagrams constitute all the possible modeled cortical responses to inhibition 
that are physiologically meaningful in Liley's MFM. 

Partition of parameter space: families reaction to anesthetics and distribu- 
tions of parameters 

Having shown the properties of Fl and F2 types of diagrams, and con- 
cluded that any of the 405 plots of the preliminary sets falls into one of the 
two, we process the whole 73,454 sets in search of global, statistically rele- 
vant, correlations between family membership and parameter distributions. 
The approximation we introduce to cope with the large number of sets in 
an automated fashion is simple: Fl and F2 are those sets respectively with- 
out and with separate lines of sn, which are divided by the exchange of sn 
branches illustrated in Figs. HC-D. Hence, Figs. [3A and[4]A.-C are considered 
as Fl, Figs. 03 and HP as F2. The second swallow tail bifurcation men- 
tioned previously is the boundary between the so-defined Fl and F2. The 
algorithm we employed searches for the position and number of cusps on sn 
lines, and assigns the family type consequently. Equilibration, continuation 
and attribution of a single set to its family takes on average 9 seconds on an 
everyday desktop computer. 

The first, global difference between families we wish to show lies in their 
response to modeled anesthetic action. According to the methods explained 
in Bojak and Liley |3], it is possible to simulate the increase or decrease 
in EEG power when anesthetics are induced, and extract the ratio between 
the power of the anesthetized cortex and the cortex at rest. Herein power 
is predicted for the system at a stable equilibrium subject to fluctuations 
induced by noise added to the thalamocortical input p ee . In line with clin- 
ical practice, they adopted the minimum alveolar concentration (MAC) of 
anesthetic agent at 1 atm pressure as the reference measure for describing 
the behavior of the anesthetized cortex. One MAC is the inspired anesthetic 
concentration needed to prevent movement in 50% of people to a noxious 
(surgical) stimulus. For example, assuming that isoflurane, a common anes- 
thetics, is employed during surgery, a patient would be maintained usually 
at 0.9 to 2.2 MAC isoflurane in an oxygen-70% nitrous oxide mixture, or at 
1.3 — 3 MAC without the nitrous oxide. According to previous studies, it is 
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assumed that 1 MAC is equal to a 1.17 vol% for isoflurane, which corresponds 



to c = 0.243 mM aqueous concentration [27], |12 



The distributions of family responses to anesthetics according to power 
ratios are depicted in Fig. |6l with F2 types being in the majority in the 
whole batch, since they turn out to be about one and half times as many 
as Fl. The histograms look similar in shape but are shifted with respect 
to one another. They are characterized by similar averages, with a relative 
difference of only about 5%, but are strikingly different at the respective 
tails. At low ratios, Fl sets have a tendency of appearing more frequently 
than F2 sets, and vice versa for high ratios. For example, for power ratio 
values smaller than 0.7, the cumulative probability for Fl is almost three 
times that of F2, whereas for ratios bigger than 1.0 the situation is reversed, 
with F2 having a cumulative probability more than three times that of Fl. 
This shows that global, dynamical patterns have an interesting degree of 
correlation with total EEG power when responding to modeled isoflurane, 
whose mode of action is representative of most anesthetics. 

Intuitively, anesthetized cortices should show a decrease in total EEG 
power when anesthetics are induced, as behavior is impaired as the cortex 
"slows down" , with power ratios at 1 MAC expected to be smaller than one. 
This is not always true, and it is not uncommon to rather observe a transient 
surge in EEG power that is clinically known as the biphasic response jij • For 
the generated 73,454 sets, a power increase larger than 1.4 at 1 MAC, has 
been taken as indicative of a biphasic response in j^], i.e., the simulated total 
power integrated over all frequencies was at least 1.4 times larger with than 
without the presence of 1 MAC isoflurane. Only 86 sets out of the whole 
batch show a biphasic response, of which around 90% are of F2 type. This is 
interesting, but the statistics are too poor to draw a general conclusion about 
the relation of F2 to the appearance of a biphasic response. It is possible that 
the predominance of F2 for high ratios is related to the cortical mechanism 
eliciting a surge in power when anesthesia is present. 

Besides responses to anesthetics, our family partition also shows that 
some frequency distributions of parameters exhibit significant differences for 
Fl and F2. In order to objectively assess these differences, we use the fol- 
lowing procedure: first, we construct parameter frequency histograms with 
n = 20 bins within the limits provided for each parameter by Tab. [TJ sep- 
arately for the two families. Two particular cases arise in this procedure. 
On one hand for and h^, which have flexible upper limits, we choose 
the maximum possible one (—65 mV) as bin limit. On the other hand, 
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occupies almost exclusively the lower part of the allowed parameter range. 
For the sake of better accuracy with the same number of bins, we set 120/s 
as the upper bin limit for 7^. Parameter sets with larger values, i.e. 12 
out of 29,131 for Fl and 24 out of 44,323 for F2, are counted in the highest 
bin. Second, we now compute the square root of the information radius (also 
known as the Jensen-Shannon divergence) between the histograms for each 
parameter j - g\ for Fl and h\ for F2 - with i — 1, . . . , n bins 



\ 



~ (ol log 2 g\ + h\ log 2 hi) - ai log 2 
with ai = - (gl + hi) , 



(10) 



which is a proper metric for the similarity of discrete probability distribu- 



tions [llj. Note that for empty bins, one defines 01og 2 = 0. Identical 



distributions have diR = 0, whereas maximal dissimilarity is indicated by 
g?ir, = 1. An example for maximal dissimilarity would be histograms g\ = 5i 



ik 



and hi = 5u with Kronecker deltas and k ^ I. 

The eight largest dissimilarities between the frequency distributions of 
the families are found for the parameters r e (0.42), a e (0.33), (0.16), 
7 ei (0.15), j u (0.15), T ie (0.12), Tu (0.11) and N? e (0.11), where the number 
in brackets is the corresponding d{ R value. These distributions are shown in 
Figs. [7] and El In the case of r e , see Fig. [7A, Fl reaches a maximum around 
0.03 s, with a slow decline for higher values, while F2 has a steadily increasing 
trend, perhaps saturating above 0.1 s. This is also the case for F2 sets in 
a e , see Fig. 03, whereas Fl gives rise to a broad unimodal distribution with 
a maximum around 4.2 mV. For both r e and cr e , the difference in frequency 
for Fl and F2 at the edges of the intervals is striking. This strongly suggests 
that r e and a e play a relevant role in selecting the global dynamics of the 
model, with a clear bias for one type of family over the other. 

Although differences in the distributions of other parameters are still 
significant, they do have a minor effect in predicting family membership, 
since the information radius distance of r e and a e is by far the largest of 
all the distributions. For example, in Fig. [7p and 7 e j in Fig. 03 retain 
separated maxima for the occurrence of Fl and F2, showing a tendency 
for larger likelihood of Fl at low physiological values and a faster decrease 
in frequency as the parameters increase. Notice how 7^ in Fig. [HA shows 
bell-shaped distributions concentrated in the subinterval 10 — 100/s, in a 
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similar configuration to that of total power ratios discussed above in Fig. El 
In contrast, for PSP amplitudes F ie , Fa, and connection strengths N?, see 
Fig. [8f3-D, we observe larger changes for F2, while Fl is closer to constant 
throughout the allowed parameter range. In these three cases at the edges of 
the physiological parameter interval one finds clear bias for one family over 
the other, but overall the difference is not as strong as for r e and a e . The 
distributions of the remaining parameters do not show marked differences 
between Fl and F2, with d{ R values between 0.012 and 0.10. 

Different oscillatory responses of families to altered inhibition 

Belonging to Fl or F2 only weakly predicts the presence of a surge or 
decrease in power when isoflurane action is simulated, since only the tails in 
Fig. [6] have a predominant presence of one family over the other. Neverthe- 
less, the variabilities of patterns belonging to Fl and F2 embody different 
responses to changes in inhibition. Thus we focus now on the oscillatory 
activity associated with Fl and F2, when changes in inhibition are modeled 
solely via alteration of the parameter R, to see how alterations in inhibitory 
PSP modify the oscillatory "landscape" families can produce. To this end 
we inspected recurrent lpar plots at fixed k, within a physiological range of 
0.75 < k < 1.25, and characterize the oscillatory activity. 

Examples of this are depicted in Figs. [9] and [10J obtained from selected 
sets in the initial batch of 405. We want to stress that the variations within 
each family are extensive: the plots proposed are only schematic indications 
and are not meant to exhaustively describe the behaviors of Fl or F2 types. 
Nonetheless, some general conclusions can be drawn, based on the most 
commonly appearing dynamical traits in the continuations in R. Firstly, 
Fl seems to be less prone to stable oscillatory activity in the physiological 
range than F2. Changes of R, and thus of T ie and Tu, tend to trigger unstable 
periodic orbits with a limited parameter span for Fl, whereas F2 shows stable 
cycles in a larger parameter range. In particular, plots where the size of the 
parameter interval of stable orbits is either very limited (Fig. [9JA.) or zero 
(Fig. [HB) are frequently observed in Fl types, whereas plots with significant 
oscillatory regimes (Fig. [HP and D) are infrequent. 

Secondly, the presence of two or more stable regimes is rarer in sets be- 
longing to Fl, as compared to F2, and furthermore they typically extend 
over more limited R ranges within the physiological interval. Hence the oc- 
currence of separate periodic orbits at the same value of R is much rarer 
in Fl than in F2. In other words, multistability is distinctive of F2. For 
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instance, the plots in Fig. [HP and D have a very limited extension of stable 
orbits for physiological values of R compared with Figure HOB. C and D. 
Since a relation between multistability and memory has been hypothesized, 
this could suggest that F2 types possess dynamical structure that facilitates 
the formation of memory. In general and on average, it appears that Fl is 
associated with a more restricted dynamical repertoire than that of F2, and 
that F2 responds to the changes in inhibition with more activity than Fl. 

The distribution of the Hopf bifurcations in different families gives some 
insight into this. In Fl, the curve of Hopf bifurcations can simply be closed, 
with no interaction with the saddle-node bifurcation. This is responsible for 
the dynamics in Fig. [UA and B. In this case, there exists a time-periodic 
solution but it is stable only on a small domain. A slightly more complicated 
situation arises if the Hopf bifurcation intersects with the saddle-node bifur- 
cation at the fh points, as in Fig. [3]A This leads to a more folded branch 
of periodic solutions, as shown in Fig. |9p and D. Again, stable oscillatory 
behavior is observed only in a small domain. In F2 instead the hb curve is 
stretched out and interacts with both branches of sn curves and generally 
with both sides of the A-shaped branch, as in Fig. [SB. This last interaction 
is, in particular, closer to the physiological range than the fh points usually 
observed for Fl. Hence, the configuration in F2 does have a greater potential 
for sustained oscillations and the coexistence of stable periodic solutions, as 
shown in Fig. HOB. 

5. Discussion 

We have illustrated how bifurcation diagrams of Liley's MFM in two pa- 
rameters, R and k, account for core effects of inhibition over the cortex and 
can be classified into two topologically distinct families. These parameters 
modify inhibitory PSPs and inhibitory self-coupling, respectively, and the 
two archetypal families are identifiable by invariant dynamical features of 
model cortex. Across a large number of parameter sets, we have listed and 
identified all the possible diagrams. The relatively small number of bifur- 
cations and combinations between branches we have found indicates that 
Liley's MFM, when it supports "realistic alpha activity" at nominal parame- 
ter values [3| , results in a relatively limited dynamical repertoire for cortical 
activity in response to parametric variations within physiologically admis- 
sible space. These dynamical structures have been classified into families 
according to their shared global topology and equivalent local properties of 
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the solutions around their bifurcations. Family membership weakly corre- 
lates with typical changes in total spectral power when the cortex is acted 
upon by anesthetic agents, but strongly with unusually large changes. At the 
same time, the likelihood for family membership is sensitive to changes in 
specific neurophysiological attributes, in particular to the standard deviation 
of the excitatory neural population firing threshold a e (i.e. the steepness of 
the neuronal firing rate function Sk) and the excitatory membrane decay time 
constant r e . As far as we are aware, these connections between dynamical 
properties of a MFM, physiological attributes of cortex and changes in EEG 
power spectra under anesthesia have never been found before. 

It is important to appreciate the generality and relative simplicity of the 
methods we have used here. We have chosen two meaningful parameters 
for the bifurcation analysis, R and k, and have then classified bifurcation 
diagrams into two families, Fl and F2, according to the absence or presence 
of separate sn lines. Family membership turned out to be separated by a 
swallow tail bifurcation. More sophisticated criteria for labeling resulting 
2par plots could be envisaged. For example, families might have been cate- 
gorized by the number, location and types of bifurcation they presented, or 
by the values of the coefficients in their normal forms. However, in our case 
the simple approach proved sufficient. The main methodological result we 
wish to underline is that by grouping the global dynamics for a representa- 
tive sample of 73,454 parameter sets of an EEG model according to agreed 
criteria, strong, distinctive and robust correlations with EEG power spectra 
on one hand and some of the physiological parameters on the other hand 
have been derived. This is novel per se, because it shows that systematic 
partitioning of complex parameter spaces according to global dynamical pat- 
terns can lead to significant results which are not otherwise accessible. This 
method seems to be one the few viable strategy for unveiling interesting re- 
lations in high dimensional nonlinear spaces and exploring typical responses 
in complex models. 

In this paper, we have focused on the action of anesthetics and ratios of 
total power, but many other developments are possible. For example, we 
could select those sets which give rise to epileptic activity or have a high 
power in the EEG gamma band, categorize their dynamical patterns and 
look for parameter correlations. Our approach is hence promising for other 
theoretical studies of the EEG, or neural activity in general. 

We want to stress that the proposed method is not dependent on Liley's 
MFM: in general, any other set of coupled, nonlinear ODEs with a compli- 
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cated but meaningful parameter space can be treated likewise. The difficul- 
ties in extracting relevant information from equations with many parameters 
and involved feedback and feedforward mechanisms should not be underes- 
timated, and our dynamical approach can offer new guidance. Here a PCA 
on the 73,453 analyzed parameter sets turns out to be completely useless. 
This is also true if we employ a PCA after dividing the sets into Fl and F2: 
the first ten principal components then account for about 55% of the total 
variance, up from 42%. Again, no valuable information on parameter rela- 
tions can be gleaned from these types of statistical methods. If we compare 
Fig. [2] with the clear separation between distributions in power ratio in Fig. [6] 
or in Fig. [7] for parameters r e and a e , we can appreciate the strength of our 
method to tease out hidden correlations. 

What is the biological meaning of family differences with respect to pa- 
rameters such as r e and cr e ? First, a small excitatory rate constant r e en- 
hances the speed with which h e (t) adapts to changes in PSP inputs, as is 
clear from Eq. |2j For r e = 0, the mean excitatory soma membrane potential 
simply becomes a weighted sum of the PSPs. Thus, low r e enslave h e (t) 
more directly to synaptic input and results in a depression of self-sustained 
oscillations, which is characteristic for Fl. a e parametrizes the local exci- 
tatory firing rate response, see Eq. [3j and hence does not directly act upon 
h e (t) dynamics, as r e does. However, a smaller <j e means that the excita- 
tory firing rate changes more rapidly for variations of h e (t), although over 
a more limited range. In the limit of a e = 0, one obtains a step function 
which instantaneously switches the population from zero to maximal firing 
at the threshold value \i e . As mentioned above, in a steady state situation 
excitatory self-feedback is essentially proportional to {N™ e + iVfJ ■ S e [h e (t)], 
and a low a e increases this self-feedback. Thus, as for r e , the reaction of h e (t) 
to changing synaptic input will be more rapid for low a e , limiting the pres- 
ence of self-sustained oscillations and hence favoring membership in Fl. One 
can then view the somewhat lesser differentiation of a e with regards to the 
families due to its boosting excitatory reactivity more indirectly. Overall, en- 
hancements in the reactivity of excitatory populations result in a diminished 
dynamical repertoire. 

These considerations on the primary role of a e and r e in influencing family 
membership are not weakened by the presence of complex feedbacks between 
neuronal populations in Liley's model, c.f. Fig. [TJ Sets of Fl types appearing 
at low r e or er e show homogeneous distributions for all the other parameters, 
with no formation of clusters or any other form of bias. This suggests that 
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the model cortex generally can be pushed towards Fl by appropriately re- 
ducing the excitatory membrane decay constant or the excitatory standard 
deviation of the firing threshold, and towards F2 by increasing them. Thus, 
we can also identify excitatory reactivity as an endogenous counterpart to 
the exogenous control. In fact, it has been shown in detail here also how 
exogenous activity, namely the excitatory thalamic input to inhibitory (exci- 
tatory) cortical populations described by p ei (p ee ), can cause changes across 
and within families. There exists a kind of mirror symmetry between p ee and 
p e i: at the increase of p e i, or decrease of p ee , we have observed a transition 
from F2 to Fl, see Fig. HJ corresponding to a transition from complex to 
simple dynamics. This indicates that the exogenous driving of inhibition (or 
a reduction in the driving of excitation) moves the cortex to dynamics which, 
on average, appear simpler. 

Thalamic input appears in Liley's MFM as an active source of control on 
the complexity of brain dynamics, having the role of modulator in cortical 
interactions. This is quite different to the classical view of thalamus as being 
the gateway through which peripherally derived sensory information reaches 



cortex[40[. For example, retinal axons project to the lateral geniculate nu- 
cleus of the thalamus, which then projects to visual cortex. In this way the 
thalamus is seen as the conduit through which all information about lower 
levels of the nervous system reaches the cortex for "processing" . Viewed from 
the perspective of our model, this would imply that thalamic input simply 
drives the cortex. However, based on our metabifurcation analysis we predict 
that such input also configures (i.e. "modulates") local and global dynamical 
responses from the cortex. 

It is now known that all areas of cortex, not just the primary sensory 
cortices (somatosensory, visual, auditory) receive thalamic input, and that 



most areas of cortex, in some reciprocal manner, innervate thalamus [35] . As a 



consequence, contemporary approaches to understanding the role of thalamus 
are now aimed towards better articulating the functional inter-relationships 
that exist between cortex and thalamus. One popular approach has been 
the modeling of integrated dynamics of thalamo-cortical activity. On this 
basis, it has been speculated that the human alpha rhythm emerges as a 
consequence of reverberant activity between cortex and thalamus. While we 
have not here considered such thalamo-cortical feedback, doing so in the con- 
text of our results leads to a number of interesting speculations regarding the 
auto-regulation of cortical dynamics: cortical feedback through the thalamus 
could initiate a sequence of transitions between bifurcation pattern families, 
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providing a way of reconfiguring cortical dynamics "on the fly" and thus on a 
time scale quite different to that associated with activity dependent synaptic 
plasticity. Unfortunately, not enough is known experimentally about thala- 
mocortical feedback at present to turn such speculations into strong model 
predictions. 

On the other hand, thalamus also appears to stimulate internal variations 
within families, since a lower value of p e i (higher value of p ee ) corresponds to 
more complex bifurcations as depicted in Fig. |5j As said, this type of 2par 
plot with two bt points has been shown to produce chaotic activity, meaning 
that "cortical state" F2 is capable of producing highly nonlinear events, and 
an exogenous enhancement of cortical inhibition tends to reduce this abil- 
ity. This and the fact that F2 types have shown to occur with the largest 
likelihood in parameter space (i.e., their total number is one and a half time 
more than Fl), stimulates some speculations. F2 may be considered as the 
cortical default state of Liley's MFM, associated with the richest dynamical 
responses to inhibition and acting as a sort of cognitive readiness state. Both 
endogenous (cortical) enhancements of excitatory reactivity and exogenous, 
subcortical increases (decreases) of inhibitory (excitatory) activity reduce the 
cortical state to much simpler behavior. 

Finally, we have implemented a partition of the parameter space accord- 
ing to the bifurcation plots. As mentioned, types of normal forms appear to 
be constant within each family and continuation lines of the diagrams are 
similar, whilst the extensions of these patterns in the Rk-space vary wildly. 
This suggest the hypothesis that a single, high dimensional master bifurca- 
tion diagram exists and governs all sets, and that changes in parameters only 
change its orientation, position and scale in the (hyper)volume. By continu- 
ing in R and k, we intersect the master diagram with a plane. The position 
of this plane determines whether we see Fl or F2 diagrams, and the two can 
be continuously deformed into each other. 

We would like to conclude with some future directions. First, an inventory 
of all the possible types of oscillatory dynamics associated with Fl and F2 
is surely attractive. How do many of them support complex events such as 
chaos and how is this reflected in parameter space? Do chaos and complex 
dynamics occur for limited, selected intervals in the parameters? The great 
hurdle is represented by the computational demand: it is still unfeasible 
to process all sets one by one and see what types of orbits are present. 
Second, a systematic study of the unfolding of the bifurcations in 2par plots 
that characterize families should be attempted. This has the potential to 
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increase our knowledge about the type of dynamics these parameter sets 
can produce. Finally, there is an interesting theoretical question concerning 
the optimization of the method we presented in this paper. What is the 
best choice of continuation parameters for achieving the highest degree of 
separation among bifurcation plots? We wonder if it would be possible to 
choose parameters or combinations of parameters such that their degree of 
clustering with respect to families is maximal. Some canonical discriminant 
analysis on the sets has been preliminarily performed, but results are not 
particularly illuminating. Achieving optimal partitioning has the potential 
to uncover biologically relevant dependencies among parameters that are still 
unknown, and to shed more light on clinically or pharmacologically relevant 
questions. 
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Table 1: Physiological ranges for spatially averaged model parameters. Common ranges does not imply shared values 
for different types k — e, i of neuronal target populations, hence physiological parameter space is 32-dimensional. For bulk 
oscillations and for absent inhibitory extracortical input, vA and r^iV^ occur only in combination: this makes parameter 
space effectively 29-dimensional. ^The upper limit for should be h r e — 5 mV, but was erroneously set to h\ — 5 mV in 
Ref. Q. Since power spectra were obtained around a chosen fixed point with h* e > for all parameter sets, this has no 
significant consequences. ^Equivalently time to peak PSP amplitude given by 1/7/&. 
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Figure 2: Principal Component Analysis of the model parameter space. Fraction 
of the total variance of parameters explained by the first ten principal components for all 
73,454 sets. Parameter space has 32 physiological dimensions, of which 29 are effectively 
distinguished here, c.f. Tab. [TJ Orange bars indicate individual contributions, brown ones 
their cumulative sum. The evident absence of strong linear correlations within parameters 
means parameter space cannot be reduced easily to a lower number of effective degrees of 
freedom. 
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Figure 3: Typical bifurcation structure of families. Continuations for two distinct 
parameter sets belonging to Fl (A) and F2 (B), respectively. Note in B that fold-Hopf 
points fh2 and fri3 of F2 occur close to the cusps C2 and Ci, respectively. The two saddle- 
node lines are shown in blue sni and green sri2, respectively, and the Hopf line hb in red. 
The number of equilibria is 3 for areas bounded by the saddle-node lines sn, and 1 outside. 
The rectangle shows the physiologically relevant range of i?fc-space. 
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Figure 4: Exogenous effects: inter-familiar transition. Metamorphosis from Fl to 
F2 through two swallow tail bifurcations, induced by decreasing thalamic excitation to 
inhibitory neurons: A Fl at p e i — 26,000/s, B Fl at p e i — 14,000/s, C just before 
bifurcation at p e i = 7,000/s, and D just after bifurcation p ej ; = 6, 800/s. The separation 
between the cusps Ci and C2 further increases for even lower values of p e i, and a pattern 
identical to the typical F2 in Fig. [3j3 appears at p e i « 6, 000/s. Analogous transitions 
occur for increasing thalamic excitation to excitatory neurons (p ee )- The overwhelming 
majority of the inspected sample of 405 sets metamorphose similarly. 
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Figure 5: Exogenous effects: annihilation of bt points. Increasing thalamic excita- 
tion to inhibitory neurons results in the annihilation of the bt points: A p e i = 3, 720/s, B 
Pei = 3, 724/s, C p e i — 3, 729/s and D p e i — 3, 735/s. The inset in A shows the triangular 
region of F2 where the bifurcations take place. Further increase of p e i shrinks hbi up to the 
coalescence and disappearance of the bt points. Note that the transformations induced by 
the thalamic input occur over a very small range of p e i and within even smaller intervals 
for R and k, demonstrating the sensitivity of the parameters. To improve visualization, 
continuation lines have a larger linewidth and the sn branch is dashed. 
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Power Ratio at 1 MAC 



Figure 6: Distributions of power ratios sorted by families. The normalized distri- 
butions of power rations for the analyzed 73,454 sets, with Fl in blue and F2 in red, and 
magnifications of their tails in the insets. There are 29,131 sets of Fl type (39.7% of the 
total) and 44,323 sets of F2 type (60.3%). The total number of bins is 250 for power ratios 
— 4, of which only — 2 is shown. Average and standard deviation of the power ratios are 
0.77 ± 0.09 for Fl and 0.82 ± 0.09 for F2. Cumulative probabilities for the distributions 
at the tails are P(power ratio < 0.7|P1) = 0.161, P(power ratio < 0.7|P2) = 0.060, and 
P(power ratio > 1.0|P1) = 0.012, P(power ratio > 1.0|P2) = 0.037. 
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Figure 7: Normalized frequencies of parameter values by families, part 1. Blue 
bars are for Fl, red ones for F2. The total number of bins is 20, within the intervals of 
physiological validity indicated in Tab.[T] Histograms A and B represent the distributions 
with the most striking differences. 
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Figure 8: Normalized frequencies of parameters values by families, part 2. See 
caption of Fig. [7] for details, ja is here binned only in the interval — 120/s, since larger 
values are extremely rare. Also note the different y-axis scales. 
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Figure 9: Continuations in R of four different Fl parameter sets at fixed values 

of k. The interval 0.75 < R < 2.00 is considered physiologically meaningful, as are the 
chosen values of k (A k = 1.25, B k = 1.25, C k = 1.00 and D k = 1.00). Green lines 
denote equilibria, blue lines correspond to stable and red to unstable periodic orbits, i.e., 
self-sustained oscillations. To avoid clutter, only the maximal h e of orbits are plotted and 
changes to the stability of equilibria between the sns and hbs are not shown. Examples of 
periodic orbits from stable and unstable branches in physiologically meaningful intervals 
are also depicted as insets: in C an unstable oscillation converging to equilibrium in red, 
and a stable one in alpha frequency range (8 — 13 Hz) in blue; in D two stable orbits 
with different amplitudes, the top one in beta (13 — 30 Hz) and the bottom one in alpha 
frequency range, respectively. Note that the latter two orbits partly overlap in R, and 
thus constitute an example of multistability. 
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Figure 10: Continuations in R of four different F2 parameter sets at fixed values 

of k. As in Fig. [3 but for F2 and with the following values of k: A k = 1.00, B k = 1.25, 
C k = 1.25 and D k = 1.25. Note how the extension of stable orbits and the existence 
of multistable regimes (i.e., multiple, different orbits occurring at the same R) is more 
pronounced compared to Fl in Fig. [9l In B and C stable periodic orbits at physiologically 
meaningful values are also shown as insets. In B, where a multistable regime is present, 
frequencies are in the beta and alpha bands for the higher and lower inset, respectively, 
whereas in C a typical alpha oscillation is depicted. The inset in D is a magnification 
showing more clearly the presence of multistability in the non physiological R range. 
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